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Abstract 

Current biological knowledge supports the existence of a secondary group of cancer cells within the 
body of the tumour that exhibits stem cell-like properties. These cells are termed cancer stem cells, 
and as opposed to the more usual differentiated cancer cells, they exhibit higher motility, they are 
more resilient to therapy, and are able to metastasize to secondary locations within the organism and 
produce new tumours. The origin of the cancer stem cells is not completely clear; they seem to stem 
from the differentiated cancer cells via a transition process related to the epithelial-mesenchymal 
transition that can also be found in normal tissue. 

In the current work we model and numerically study the transition between these two types of 
cancer cells, and the resulting “ensemble” invasion of the extracellular matrix. This leads to the 
derivation and numerical simulation of two systems: an algebraic-elliptic system for the transition 
and an advection-reaction-diffusion system of Keller-Segel taxis type for the invasion. 

Cancer research aims to understand the causes of cancer and to develop strategies for its diagnosis and 
treatment. In this effort, mathematics contributes with the modelling of the biological processes, and the 
corresponding analysis and numerical simulations. As a research field, it has been very active since the 
1950s; see e.g. [23, 3, 8] and has spanned over a wide range of applications from intracellular bio-chemical 
reactions to cancer growth and metastasis, see e.g. [25, 24, 2, 11, 9, 21, 30, 22, 13, 32, 12], 

Intrinsically, cancer cells exhibit higher proliferation rates than normal cells. Recent evidence though 
have revealed the existence of a subpopulation of cancer cells that posses lower proliferation rates than 
the rest of cancer cells, exhibits stem-like properties such as self-renewal and cell differentiation, and 
are able to metastasize i.e. to detach from the primary tumour, invade the vascular system, and afflict 
secondary sites [4, 19]. These are termed Cancer Stem Cells (CSCs) and constitute a small part of 
the tumour, while the bulk of the tumour consists of the more proliferative Differentiated Cancer Cells 
(DCCs), [31, 27] 

The origin of CSCs is not completely known; a current theory states that they are de-differentiated cancer 
cells originating from the more usual DCCs [12, 27]. This type of transition of cancer cells, influences 
their cellular potency and seems to be related to a type of cellular differentiation program that can be 
found also in normal tissue, the Epithelial-Mesenchymal Transition (EMT) [19]. The EMT takes place 
in the first step of metastasis, the invasion of the Extracellular Matrix (ECM) by the primary tumour, 
cf. [31, 14], 

The EMT is triggered by the micro-environment of the cell [26]. The protein, for example, Epidermal 
Growth Factor (EGF) and its corresponding cellular receptors EGFR plays pivotal role in this triggering 
[33, 28]. This process can be briefly described as follows: the EGFs function as “on”/“off” switches 
that stimulate cellular growth and proliferation. In certain pathological mutations, they get stuck in the 
“on” position and cause unregulated cell growth, untimely EMT, and over-expression of EGFRs. The 
over-expression though, of EGFRs has been connected to the tumor formation in mammary epithelial 
cells, see e.g. [10]. 
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From mathematical perspective, the modelling of ECM invasion takes into account a special type of 
taxis, namely haptotaxis. In haptotaxis, the cells move biased, along directions being dictated by the 
availability or the gradients of either extracellular adhesion sites (which are located on the ECM), or of 
ECM-bound chemoattractants. 

From the macroscopic point of view haptotaxis (as the other forms of taxis) is typically modelled under 
the paradigm of Keller-Segel systems which was first introduced in [15] to describe the aggregations of 
slime mold. Since their first derivation, they have been successfully applied in a wide range of biological 
phenomena spanning from bacterial aggregations to wound healing, and cancer modeling. In cancer 
growth modelling in particular, Keller-Segel type systems were extended to include also enzymatic in¬ 
teractions/reactions, yielding Advection-Reaction-Diffusion (ARD) models, see for instance the generic 
haptotaxis system proposed in [2] (promoting ECM degradation by Matrix Metalloproteinases (MMPs)), 
or the chemo-hapto-taxis system proposed in [5] that promotes the role of the enzyme plasmin. 

From a numerical point of view, the treatment of this type of problems is a challenging task, since their 
solutions develop (generically) heterogeneous spatio-temporal dynamics in the form of merging/emerging 
clusters of cancer cell concentrations, see for instance Figure 3. To treat these dynamics we employ 
a second order Finite Volume method equipped combined with a third order Implicit-Explicit time 
integration scheme, see Section 2 and our previous works [17, 16] for more details. 

In this work, our objective is twofold; first we model the EMT transition from DCCs to CSCs, triggered 
by the availability of free EGF molecules in the extracellular environment, and of the corresponding 
EGFR cellular receptors, and secondly we embed this EMT model in an ECM invasion system that 
addresses the “ensemble” CSCs and DCCs. We subsequently perform some numerical experiments to 
investigate properties of the combined system. 

In more detail, we construct a parabolic system describing the concentrations of EGFs and EGFRs 
and their attachment/detachment dynamics. These processes are assumed to be much faster than the 
cellular motility and invasion of the ECM, hence we rescale the previously developed system and deduce 
an elliptic system that addresses the free and bound EGF. The rescaled elliptic system serves as the 
large time asymptotic limit of the parabolic system. 

To model the ECM invasion, we employ an ARD haptotaxis model that we adjust to include the CSC- 
DCC “ensemble” and the EMT between them. In particular we assume that the cancer cells diffuse in 
the environment, move responsively to the gradients of the ECM, proliferate in a logistic manner, and 
most notably undergo EMT. We also include the ECM degradation by a matrix degenerating MMP, 
which in turn is assumed to be produced by both families of cancer cells. 

With the combination of these models we are able to reproduce the current biomedical understanding 
that the CSCs a) stem from the bulk of DCCs, b) they separate themselves from the main body of the 
tumour, c) they invade the ECM while exhibiting highly dynamic clustering phenomena. 

Similar ECM invasion systems that can be found in the literature, take into account multiple subpopula¬ 
tions of cancer, can be found (see e.g. [1, 7]). These systems however control the transition between the 
subpopulations in terms of time dependent Heavyside step functions, while in our approach the suggested 
models connect the microscopic dynamics of EGF with the macroscopic dynamics of EMT and ECM 
invasion. 


1 The EMT coefficient 


Typical modeling of the ECM invasion by a single type of cancer cell includes terms for diffusion, taxis, 
reaction, and proliferation. For two types of cancer cells, the corresponding equations apply with an 
addition of a coupling term, which in our case describes the transition from DCC to CSC via EMT, i.e. 


dc CSG 

—-— = diffusion + taxis + uemt c DCC + reaction + proliferation 
at 

dc DCC 

——— = diffusion + taxis — /temt c DCC + reaction + proliferation 


(1) 
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where c csc , c DCC represent the densities of the CSCs and the DCCs respectively. The coefficient memt 
describing the probability of the EMT transition is given by Michaelis-Menten kinetics (see also [34]) 

[ff] DCC 

Memt = Mo- , r iDCC ’ ( 2 ) 

M 1/2 + 

where [< 7 ]P CC represents the concentration of the EGF molecules that are bound on the corresponding 
EGFR receptors on the membrane of the DCC cells, and where Mo is the maximum EMT rate, and M 1/2 
is the Micliaelis constant, i.e. the bound EGF needed to generate half maximal EMT. 

In addition to [g]]C cc , we set \d)f to represent the concentration of free EGF molecules (in the extra¬ 
cellular environment), [<?]o the total (free and bound on both DCCs and CSCs) concentration of EGF. 
Accordingly, [r] 0 , [r]b, and [r]/ denote the total, occupied, and free receptors on the surface of both types 
of cancer cells. 

Throughout this work we assume that all [g]_ and [r]_ concentrations depend on time and on space; 
this is so since the free EGF molecules are able to translocate (for example diffuse) in the extracellular 
environment, and the bound EGFs and the EGFRs are to be found on the surfaces of the (also time and 
space dependent) CSCs and DCCs. We assume moreover that the total EGFs and EGFRs obey local 
mass conservation properties, i.e. 


( 3 ) 

( 4 ) 

for every t > 0 and ifSlcl bounded. 


[ 9 ]o(t,x) = [g]f(t,x) + [g] b (t,x) 
[' r]o(t,x ) = [r\ f (t,x) + [r] b (t,x) 


The parabolic system. The dynamics dictate that free EGF molecules bind onto free EGFR receptors 
with rate k+, while at the same time bound EGF molecules detach with rate k- (see also [34]). It is 
assumed that the k+ and fc_ rates are the same for both types of cancer cells. That is, the time evolution 
of both the free and the bound EGF molecules is given by the system: 

I d T [g]b= k+[g\f[r\ f - k-[g] b , 

\ 9r[a}f = D f A [g]f- k+[g]f[r\f + k-\g\ b . 

where r is the corresponding time variable. We have also assumed that the free EGF molecules diffuse 
in the extracellular environment. 

Recalling now the mass conservation of the EGFRs (4), we can deduce that: 

M/ = Wo - Mb 

= A csc c csc + A DCC c DCC - \g} b , (6) 

where A csc , A DCC represent the total “number” of EGFR receptors per cell per family, and where we 
have used that the bound EGF molecules coincide with the bound EGFR receptors. 


The elliptic system. The EGF dynamics (i.e attachment/detachment onto EGFRs and diffusion) are 
much faster than the dynamics of EMT and ECM invasion. We rescale hence the time variable r in 
System (5), to match the time scale t of the systems (1), (19). To this end, we first rewrite (5) in the 
form 

d T u = DAu +p(u), 

where u = ([<?]h, [g\f) T , D = (0, Df) T , and where p denotes the nonlinear reaction terms in (5). We 
employ now the change of variables 

v(t,x) =u(T(t),x), (7) 

with the “much faster” time variable r given by 

r{t) = ~, 0 < e « 1 . ( 8 ) 

£ 
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( 9 ) 


Relation (7) yields after differentiation as 


ed t Y = D Av + p(v), 


since 

Ay (t, x) = Au (r{t),x), p(v(t, x)) = p(u(r(t), x)), ^ = 0. 

at e at 

We denote the above introduced variables in the t timescale as follows 


f 9 b(t,x)\ 
\g f (t,x)J 


f[g] b (T(t),x)\ 


= Y (t,x). 


Taking now in (9) the formal limit as e —> 0, we deduce the elliptic system that corresponds to (5): 

f 0 = k+g f ( A csc c csc + A DCC c DCC - g b ) - k.g b , 

\ 0 = D f Ag f - k + g f (X GSG c csc + A DCC c DCC - g b ) + k.g b . 


or 


0 = ± gf{ X GSG c GSG + a dcc c dcc - g b ) - g b , 
0 = A g f . 


( 10 ) 


where ko = — represents the attachment/detachment ratio of the EGFs onto the EGFRs. Though (10) 
does not include time derivatives, both concentrations g b and gf depend on time t through the densities 
of the cancer cells c CSG and c DCC . 


Neumann and additional conditions We augment the system (10) with homogeneous Neumann 
boundary conditions, 

%M = 0 , (ii) 

on 

Remark 1. The system (10)-(11) has no unique solution. This can be seen by making the ansatz 
g b = K(X csc c csc + X DCC c DGC ), K G (0,1). For this choice the constant function 

, K 

0! = ^TTk- 

satisfies the algebraic equation of (10), the boundary condition (11), and Agf = 0. Thus for every K, a 
classical solution of (10), (11) exists. 

To recover uniqueness in the system (10)-(11) we impose an additional condition, namely we set the total 
mass of EGF to be constant (in time t), 

M= g 0 (t,x)dx= / g b (t,x)dx+ / gf(t,x)dx. (12) 

J ■/ Q o O 

Condition (12) comes as a natural consequence of the system (5) endowed with zero Neumann boundary 
conditions. The total mass of EGF in this system does not change with time r since, after integration 
of the second equation of (5) over f2, we see that 

- 7 - [\g]bdx + ^- f[g\ f dx= I d T [g\ b dx+ [ —d T \g\ b + D f A[g\ f dx = 0. 

Jn d T Jci Jn 

In effect, the total mass of EGF in the elliptic system (10) given by (12) will be equal to 

M = [ [g]° f dx+ [ [ g] b dx , 

J Q, J £1 

in the parabolic system (5) for given initial concentrations [g]® and [g]°. 
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Figure 1: Plot of pemt for the system (19)-(22) at time t = 6. Note that it is space dependent and that it 
“lives” over the support of the DCCs. 


Solution of the elliptic system. We note at first that the Laplace equation in System (10) has a 
classical solution for which the gradient Vg/ is unique. We can thus deduce that gf = gf is constant. 
Accordingly, the first equation of the system (10) implies 


9b 


kp 9f ^CSC C CSC ^DCC c DCC) 
1 + T -Qf V 

(A 


+ kE~9f 
9f 

kD + gf 


,CSC C CSC _|_ yDCC c DCC) 


(13) 


Similar, but slightly more involved, is the derivation for the DCC-bound EGF. We repeat the assumption 
that the attachment/detachment rates fc+, fc_ are the same for both families of cancer cells, and that the 
superposition principle [r]/ = [r]^ sc + [r]^ 00 with [r]^ sc , [r]^ 00 the free EGFRs on the CSCs, DCCs 
respectively, holds. Accordingly the DCC-bound EGF is given by 


DCC 

9b 


9f ^DCC„DCC 

kD + gj 


(14) 


On the other hand, the free EGF c/f (which is constant) satisfies the quadratic equation 

0 = |D| ±-g 2 f +(\n\ + -L [ A CSC c CSC + A DCC c DCC da . _ g _ M (15 ) 

kn J V k D J n k D ) 

obtained after integration of (10) with respect to x and taking (12) into account. The equation (15) has 
a single positive root since, 

A = (|fi| + -! [ A csc c csc + A ucc c DGG dx - + 4-^ |fi|M > 0 (16) 

\ ko Jsi ko J kD 


and 


The root itself i.e. the free EGF, is given by 


kaM 

T 


< o. 


_ (|Q| + ^ ; n A csc c csc + A DCC c DCC da . _ M _ l) + v^A 

2^ 

- (|fi| k D + f n A csc c csc + A DCC c DCC ^ - M ) + k D V A 

2|D| 


(17) 


Remark 2. We note that gf, though constant in space, is not constant in time since f n \ csc c csc + 
A DCG c DGC dx 

varies with time t. In effect, gjf >CG and hemt depend on time. Note also that although gf 
is constant in space, g GGG (14) is not constant in space, so neither is pemt (2) constant in space. 
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(a) t = 0 (b) t = 2 



(c) t = 40 


Figure 2: Numerical solutions of (5) converge towards the numerical solution of (10) as t increases, if 
Neumann boundary conditions (11) are assumed for both systems and if system (10) is augmented by the 
condition f n gb+gf dx = Jq[s]/ + [<7]6 dx. The parameters Df = 0.4, fc+ = 0.5, k- = 0.3 and time independent 
x csc c csc + x dcc c dcc = e -io(,-i) 2 for x 6 n = (0 , 5 ) were used. 


Combining now (2) with (14) we deduce that 


g f A DCC c DCC 

MEMT 1^0 j ~ ~ . nnn DOC 1 ’ 

A<i/2 k D + 1*1/29 f + g/X u ^c u ^ 


(18) 


which is non constant in space. We refer to Figure 1 for a graphical representation of /temt hr the 
experimental scenario (19)-(22) at time t = 6. 


Comparison of the systems (5) and (10). The benefit that stems from the lack of time variable 
in (10) is the instantaneous “sensing” of the cancer cell concentration as opposed to the “large” time 
asymptotic convergence of the system (10). 

The solutions of the elliptic system (10), coincide with the large time asymptotic solution of the parabolic 
system (5) endowed with zero Neumann boundary conditions (11) if a selection criterion in particular a 
prescribed total mass of EGF M = / n [g]° dx + / n [sr]° dx in (10)- is employed. This can be also observed 
numerically in Figure 2. 

Hence, at every instance of the “large” time scale t, at which EMT and cell movement take place, we 
can employ the elliptic system (10) to compute the density of the free EGF (17) and accordingly deduce 
the EMT coefficient from (18). 

From a numerical point of view, this process is cheap since no PDE solvers are needed. 


2 Invasion system 

The ECM invasion model we propose is based on the “generic” invasion model of [2], In this respect 
it includes haptotaxis, degradation of the ECM by MMPs, and production of the MMPs by the cancer 
cells. In addition to this basic model, we assume a secondary family of cancer cells, EMT transition 
between DCCs and CSCs, logistic proliferation for the cancer cells, and production of the MMPs by 
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(a) t = 0.5 


(b) t = 1.6 


(c) t = 6.0 




Figure 3: Time evolution of the cancer invasion system (19)-(22). (a) the CSCs are produced -via EMT- 
by the DCCs, (b)-(c) due to their higher motility, the CSCs escape the main body of the tumour and invade 
the ECM faster than the DCCs. Clearly two propagation fronts have been formed, (d)-(f) the CSCs present 
dynamic merging and emerging concentrations; the DCCs propagate into the ECM slower, smoother, and 
without the wealth of phenomena of the CSC. 


both families of cancer cells. It reads itself as follows: 


dc CSG 

dt 

dc DCC 

dt 


= DrVc csc - XiV (c GSC Vu) + memtc dcc + W c CSC (l - c csc ) 
= D 2 Vc DCC - X2 V (c DCC Vn) - M emtc dcc + ^ cc (l - c DCC ) 


dv r 
— = —omv 
at 


dm 

dt 


= D m Am + aic csc + a 2 c DCG — /3m 


For the rest of this work we assume that (19) is equipped with the initial conditions: 

'c o csc (a;)=0, 

= e — ^05 

0 [) ’ , xed, 

v(0,x) = l-0.5co CC (a;), 

m( 0,x) = 0.2 c® GC (x) 

where 12 = [0, 7.5]. The parameters employed are: 

V = {D 1 = 3.8 x 10“ 4 , D 2 = 3.5 x 10“ 4 , D m = 2.5 x 10“ 3 , Xi = 4 x 10~ x , x 2 = 4 x 10" 2 , 

fi\ = 0.2, fj, 2 = 0.1, a\ = 0.5, a 2 = 0.5, 8 = 2, (3 = l}, 

and for computing /temt, we have used 

^ EMT = {Mo = 0.55, A csc = 1.4, A dcc = 1, M = 1, k D = 0.5, /r 1/2 = 2}, 


(19) 


( 20 ) 


( 21 ) 

( 22 ) 
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Figure 4: Propagating fronts-distance (left.), and total masses (right) of the CSCs and the DCCs as functions 
of k,D = -r— (upper), M (middle), and /ij /2 (lower) in the case of the cancer invasion system (19)-(20) with 
the rest EMT parameters as in V^. EMT , (22); the ECM parameters are as in V, (21). An increase of both, the 
ratio kD and the Michaelis constant weakens the migration and the growth of CSCs. Increasing the EMT 
mass M on the other hand yields a larger CSC/DCC distance and a higher amount of CSCs. 


see also (18), (17), (16). 

The above parameters in V were chosen in the magnitudes in which they were used in similar models (cf. 
[1], [20]). Due to the properties of cancer stem cells, we chose their diffusion, and haptotaxis coefficients 
to be larger than the corresponding numbers for the differentiated cancer cells. The parameters in (22) 
were extracted by numerical experimentation and are such that the system (19) exhibits the current 
biological understanding of the problem. 

The numerical method employed on this system is a second order Finite Volume (FV) method in space 
combined with a third order implicit-explicit time integration scheme (IMEX3). This numerical method 
has been developed for such type of cancer invasion models that exhibit merging/emerging concentration 
phenomena. For more details we refer to our previous works [17, 16] and to [6, 18]. 

In Figure 3 we present the dynamics of the system (19)-(22). After a short period of time, some DCCs 
at the ECM degenerating tumor front undergo de-differentiation and become CSCs. They exhibit higher 
motility and invasiveness than the DCCs and escape the main body of the tumor, they also exhibit more 
dynamic phenomena than the DCCs. 




































We further study the influence of modifications of the EMT coefficient on the solution of the system. 
Therefore we compute numerical solutions of the system (19)-(21) and various V mMT , and measure the 
distance (in x) from the propagating CSC front to the slower propagating DCC front at t = 20 as well 
as the total mass of CSCs and DCCs on f1 at the same time instance. 

In Figure 4 we first present the effect that various values of ku have in the afore mentioned propagating- 
front-distance and CSCs, DCCs masses. An increase of the detachment/attachment ratio yields a de¬ 
crease of the front distance and the mass of CSCs, but an increase of the total mass of DCCs. Similarly, 
Figure exhibits the effect that M has on the distances of the fronts and the masses. Both, the front 
distance and the CSC mass increase with M but the ascent stagnates for large M. The opposite is true 
for the mass of DCCs. Finally, Figure 4 shows that the model proposes less DCCs and a smaller distance 
of the fronts for larger values of [i\i 2 . 

This behavior matches our biological understanding of the EMT system. A larger amount of EGF (large 
M) strengthens EMT and thus proliferation and invasiveness of CSCs. An increase of detachment over 
attachment of EMT molecules (fco) decreases the number of bound EGF and weakens the EMT process, 
which is also the case when heightening the Micliaelis constant ii\i 2 . The saturating behaviour at high 
M , and low kp> and /X 1 / 2 , respectively, reflects that maximal receptor saturation and /iemt, respectively, 
is reached. Thus, further increase/decrease of parameters does not change the dynamics, and distance 
and proliferation of CSC any more. 


3 Conclusion 

In this work we give a mathematical insight to the EMT and its impact on the process of cancer invasion. 
In particular, we develop an algebraic-elliptic model (10) to compute the amount of free EGF molecules 
in the vicinity of a tumour. We are then able to evaluate the non-constant /xemt coefficient that drives 
the de-differentiation from DCCs to CSCs. Subsequently, we propose a Keller-Segel type system (19) for 
the ECM invasion by both types of cancer cells. 

The numerical simulations that follow, although they serve here as proof of concept, exhibit clearly that 
the “combined” system is able to reproduce the current bio-medical understanding regarding the EMT 
first step of cancer metastasis. 

Our future work includes the extension of the EMT and the ECM invasion models to 2-dimensional 
spaces, the inclusion of further bio-chemical interactions of the cancer cells with the environment, and 
the use of more biologically relevant parameters. 

The effort invested in this work, and its importance is justified by the fact that in the evolution of 
cancer, the regulation of EMT is seen as potential drug target in cancer medicine [29, 10]. A better 
understanding, hence, of these dynamic transitions and their relation to the cancer invasion of the ECM 
will be beneficial to optimize drug targeting. 
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